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In this paper, we propose a nonlinear transport model for an aviation network. The takeoff 
rate from an airport is characterized by the degree of ground congestion. Due to the effect of 
surface congestion, the performance of an airport deteriorates because of inefficient configurations 
of waiting aircraft on the ground. Using a simple transport model, we performed simulations on a 
U. S. airport network and found a global jamming transition induced by local surface congestion. 
From a physical perspective, the mechanism of the transition is studied analytically and the resulting 
aircraft distribution is discussed considering system dynamics. This study shows that the knowledge 
of the relationship between a takeoff rate and a congestion level on the ground is vital for efficient 
air traffic operations. 


I. INTRODUCTION 

Faced with rapidly growing demand in air traffic, re¬ 
searchers have addressed aggravating congestion prob¬ 
lems from various perspectives M- Airport congestion 
related to departing aircraft on the ground, i.e., surface 
congestion, is known to be a serious problem, and its rela¬ 
tionship to airport performance has been recently investi¬ 
gated @, [9]. Figure Oja) shows the empirical takeoff rate 
at Philadelphia International Airport presented by Sima- 
iakis and Balakrishnan Q. When the airport is not con¬ 
gested, increasing the number of aircraft on the ground 
directly leads to a high takeoff rate. However, after the 
number reaches a critical point, the performance dete¬ 
riorates due to inefficient aircraft configurations. Note 
that such unimodal relationships can be universally ob¬ 
served in the motion of self-propelled particles (lOl . I 111 
(e.g., vehicular traffic HI, HU, pedestrian flow |l4 ll5(| . 
and interacting particle systems muni). 

In this paper, we aim to disclose the physical aspects 
of air traffic from the surface congestion viewpoint. For 
this purpose, we consider a network of airports (nodes), 
the performance of which is determined by their conges¬ 
tion levels, as shown in Fig. mb). Note that network 
approaches to describe_air traffic have been proposed in 
previous studies [Tsi - liiot : however, these studies did not 
focus on the surface congestion effect. In addition to 
the conventional approaches in the field of management 
science, the network modeling of aviation transport has 
recently been attempted (HU, |22| from a network science 
HH-dH perspective, supported by the recent understand¬ 
ing of airport networks [27f43 lj| . However, the knowledge 
of congestion dynamics from a network perspective is still 
limited. 

This study presents a simple model for aviation trans¬ 
port networks based on the density-flux relationship lo- 
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cally defined in each airport. Our primary goal is to 
evaluate the stability of global air traffic against distur¬ 
bances. Because of assumptions imposed in the model, 
we cannot expect a quantitative description of the real 
systems. However, taking advantage of its simplicity, we 
aim to clarify the fundamental features of a spontaneous 
jam formation. 

The remainder of this paper is organized as follows. 
The next section defines the model and simulation con¬ 
ditions. In Section uni we show a spontaneous jamming 
transition phenomenon and its dynamics. Finally, in Sec¬ 
tion El we discuss the implication of these phenomena 
and future work. 



FIG. 1. (a) Empirical relationship between the takeoff rate 
and surface congestion at Philadelphia International Airport 
(data obtained from Ref. ||.). The shaded area indicates the 
standard deviation of the data, (b) Flow between two nodes 
in the network model. 


II. MODEL 

A. Problem formulation 

Consider N nodes, labeled 1, • • • ,N, representing each 
airport. The connectivity of the nodes is defined by an 
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adjacency matrix A, whose element Aj 3 indicates the 
presence (A i; j = 1) or absence {A i3 = 0) of a link be¬ 
tween the two nodes i and j (i,j = l,--- , iV). Time 
development of the number of aircraft in node i , rii , is 
described by the equation of continuity as follows: 

8 N N 

-^i = Y J A ji J{p j )~Y, A ^P^ (!) 

j =1 j =1 

where pt and J denote the density in node i and out- 
flux from the node, respectively. Here, for the sake of 
simplicity, we assumed that each link has an identical 
capacity. By assuming that each node has a capacity 
that is proportional to the degree, ki = and 

defining pi = rii/ki [321, the time evolution of a density 
set, {pi, • • • pat}, is written as 

8 1 N 

= (2) 

* j=i 

In the context of vehicular traffic, this formulation is 
known as the cell transmission model proposed by Da- 
ganzo [33l - l35| , in which each road segment defines a node. 
Previous studies focused on a simple graph comprising a 
small number of nodes, and very little is known about 
the dynamics on complex networks. Note that Sun and 
Bayen [36j] applied this cell transmission model to air 
traffic in a different manner, in which nodes correspond 
to discretized aviation routes. 


B. Simulation conditions 


We use A, j defined from the U. S. airport network [37| . 
The data comprise 500 major airports and 2980 undi¬ 
rected links (A.jj = Aji ), and the network is connected 
(i.e., there is at least a single path of links between any 
two different nodes). In this paper, we define an explicit 
form of flux-density relation as follows: 


J(Pi) 


Pi{ 1 Pi ) if Pi ^ Pmax 5 

Pmax(l Pmax) if Pi ^ Pmaxj 


( 3 ) 


where p max is the maximum density that a single air¬ 
port can accommodate on the ground. Above this den¬ 
sity, aircraft have to stay in the sky near the airport 
(they still belong to the destination node); J does not 
decrease since the congestion level on the ground remains 
the same. Note that in Fig. Oja), only aircraft on the 
ground are counted (i.e., p m ax corresponds to the right 
end of the figure.). We assumed a quadratic function 
to describe this part as the simplest smooth unimodal 
function. This function is often used to model dynamics 
in which capacity constraints play a role (e.g., vehicular 
traffic [13, interacting particle systems [ItJ , and popula¬ 
tion dynamics [38}). 

This function is illustrated in Fig. Ufa), hr which we 
set pmax = 0.7. We define the critical density, p cr , as 


the density that gives maximum J. In the function (|3|) . 
Per = 0.5. 

At t = 0, the initial density of nodes are set al¬ 
most uniformly as p, which is slightly fluctuated with 
a random variable 5 drawn from a uniform distribution 
[—0.005,0.005]. Each simulation is performed until the 
time development of densities stops. We study the char¬ 
acteristics of the stationary state by investigating key 
quantities, which will be defined subsequently. Unless 
otherwise noted, results are averaged over 1000 trials. 


III. RESULTS 


A. Global flux-density relationship 


We define the global flux, J, as the average flow per 
link as follows: 


j(p) = 


l^i =1 2^j =1 



AjApr) 


( 4 ) 


where pf 3 denotes the density pt in a stationary state. 
When the density is uniformly distributed in the station¬ 
ary state, J(p) coincides with the local flow J(p) at t = 0. 
However, as illustrated in Fig. ET a), one can find a gap 
between J(p) and J(p) for a certain range of p due to the 
nonuniform distribution of pf 3 . For p < p cr = 0.5 and 
P > Pmax = 0.7, the uniformity of pf 3 = p is conserved 
for all nodes, whereas for the middle initial density, some 
nodes take a single constant value of pi = pi while the 
others distribute in the range of pi > p ma x [Fig. [21(b)]. 
The lower bound density, p/, is the density value that 
satisfies J(p z ) = J(p m ax) (pz < Per) [see Fig. E^b)]. In¬ 
terestingly, at the critical density (p = 0.5), only a small 
number of nodes draw density from neighboring nodes. 

This jamming transition can be understood as a bifur¬ 
cation phenomenon. 


B. Global jamming transition 

Here, we investigate the stability of the system state 
that has uniform density distribution, pi = p {i = 
1, • • • ,N), by considering the behavior of a small per¬ 
turbation e; [we also define a vector, e = (ei, • ■ ■ , £n) T ] 
added to the state as pi = p + e*. Linearizing Eq. © 
yields 

8c 1 N 

~q£ = £7 J '(P) A H e i - J \p) e i’ (5) 

1 3 =1 

JU = - J'(p)Me, (6) 

where matrix M is defined by = Aji/ki {i ^ j) 
and Mu = — 1. The stability of the system is directly 
connected to the absence of the positive eigenvalues of the 
matrix, —J'(p)M. By applying the Perron-Frobenius 
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FIG. 2. (a) Local (J(p); gray) and global (J(p); red) fun¬ 
damental diagrams, (b) Stationary density distribution of a 
single trial. Each circle represents the density in a single node, 
(c) Global correlation coefficient. When the initial density p is 
small and large, all nodes have the same density value, which 
is reflected in the agreement of the local and global funda¬ 
mental diagrams [panel (a)]. On the other hand, when p is 
in the middle regime (pi < p < p max ), the uniformity of the 
density is broken [panel (b) and (c)[, and the global flux is 
deteriorated [panel (a)]. 


theorem to M -\-I (I is a unit matrix), one can find that 
its eigenvalues A* (i = 1, • - - , N) are not larger than 1. 
The eigenvalues of the matrix M are given by A — 1. Note 
that by definition, A, M, and M + I are all irreducible, 
and each element of M + J is not less than zero; thus, 
the Perron-Frobenius theorem can be applied. Therefore, 
the stability is dependent only on J'(p), i.e., the system 
is unstable when J'(p) < 0. Thus, local jamming in the 
fundamental relationship J(p) directly triggers a global 
jam. 

The flux deterioration shown in Fig. Ela) can be un¬ 
derstood as follows. When the system is not congested 
[Fig. [3][a)], pi = p is locally stable and the density is 
uniformly distributed to all nodes. Conversely, in high 
density cases [Fig. [H[b)] this state is no longer stable, 
and the flow between nodes equilibrates at a lower den¬ 
sity value pi = pi and higher density value pt > p max . 
Hereafter, we refer to this state as the globally jamming 
state. 


C. Global density correlation 


In the previous subsection, we found that the density 
distribution is imbalanced. To evaluate the regularity in 
this stationary density distribution, we define an indica¬ 



FIG. 3. Stability of (a) subcritical and (b) supercritical den¬ 
sity. 


tor variable cq, which reflects the density in node i: 

J 1 if Pi > Per, 

[-1 if Pi < Per¬ 
using this variable, correlation between two neighbor¬ 
ing nodes is quantified as , which returns 1 if two 

neighboring nodes have the same a and returns —1 other¬ 
wise. Figure [2][c) illustrates global regularity of the den¬ 
sity distribution defined as the average correlation degree 
per link 


C = 


£f=i ** 


( 8 ) 


where (x) denotes an ensemble average of x. Correspond¬ 
ing to the aforementioned bifurcation, C is constant 
(C = 1) for p < p CI and p > p m ax, and in the middle 
density regime, it suddenly decreases to approximately 0, 
explaining the breakage of the uniform state. From this 
value we can conclude that the global regular structure 
is absent in terms of density correlation in the globally 
jamming state. 


D. Local density correlation and dynamics 

Next, we focus on the local regularity of the stationary 
density distribution. Here, we set p = 0.6 and investigate 
the average local correlation per node 



As shown in Fig. HD a), c, is highly dependent on fc 7 > 
Decreasing values of k correspond to values of Ci decreas¬ 
ing to near — 1, which indicates that a node tends to 
have an opposite a against its neighboring nodes with 
high probability. However, as the number of neighboring 
nodes increases, it becomes difficult to take the opposite 
value against all the neighboring nodes, and thus the lo¬ 
cal correlation vanishes. In other words, the presence of 
high-degree nodes reduces the absolute value of the local 
and global correlations in the globally jamming state. 
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FIG. 4. (a) Local correlation coefficient, (b) Average density 
distribution. Each plot corresponds to averaged stationary 
density in one node. The inset in panel (b) shows degree 
distribution of the focal network. 

In addition, the density also depends on ki [Fig. |U[b)]. 
In hub nodes (large nodes), the average density is ap¬ 
proximately 0.6, which coincides with p, whereas in nodes 
with small ki (small nodes), the average density is sig¬ 
nificantly larger. This can be explained by the effective 
relaxation speed of the density [Figs. [SJa) and (b)]. This 
dependency is supported by Fig. [5](b), which shows the 
average of saturation time T) = min{f|/?,(f + 1) = 

This is because, for a large node, the density inflow to 
the node is averaged over a large number of neighbor¬ 
ing nodes and its speed of variation is relatively slow. 
Therefore, at the beginning of each simulation trial in 
the middle density regime, the change in density is dom¬ 
inant in small nodes, which is strongly influenced by the 
initial perturbation 5: for positive S, the density starts 
to increase, and for negative S density starts to decrease. 
Thus, one half of the small nodes become congested and 
the other half loses density. Since the drag of density is 
faster than discharge [Fig. [5](a)], and there is no upper 
limit on the density, the average final density in small 
nodes becomes large. Conversely, in large nodes, at the 
beginning, it tends to be deprived of density neighboring 
small nodes, which leads to a less probability of conges¬ 
tion occurring in these nodes. As a result, the density in 
large nodes is suppressed. 


IV. DISCUSSIONS 

We proposed a transportation model for aviation net, 
focusing on the surface congestion effect. The model ex¬ 
hibits a jamming transition, which is connected to a bi¬ 
furcation of the system equations. Interestingly, when 
each component of the network (airport) is jammed, char¬ 
acterized as a negative differential coefficient of a flux- 
density relationship, the uniform density distribution be¬ 
comes unstable, decreasing the total network flux (global 
jamming). Through relaxation from a uniform to an im¬ 
balanced distribution, densities on two neighboring nodes 
tend to vary inversely, i.e., one increases and the other 


FIG. 5. (a) Sample time series of density for large (k > 94) 
and small (k = 1) nodes, (b) Average relaxation time of each 
node. Each plot corresponds to averaged relaxation time of 
one node. 


decreases. However, on a node having large number of 
links, this correlation vanishes, breaking global correla¬ 
tion structures. In this model, density in a large node 
changes more slowly than that in a small node because 
the density inflow to the node is averaged over many 
nodes and does not change sharply. Note that because 
the model assumes that the capacity of each node is pro¬ 
portional to its degree, the time constant of the density 
variation is normalized against the degree. Because of 
the difference in time coefficient, small nodes are more 
congested on average. In a small node, however, there 
may be small traffic demand, which contradicts the as¬ 
sumption that each link has the same amount of traffic. 
Hence if we take the traffic capacity of links into consid¬ 
eration, the time constant in small nodes is reduced and 
thus the final density distribution becomes more uniform. 

On the other hand, the global jamming transition pre¬ 
sented in this study is not influenced by these simplifi¬ 
cations because it is caused by a fundamental bifurca¬ 
tion observed even in a two-component system }.35j . In 
previous studies, this type of jamming transition in air 
traffic systems has rarely received attention. An effective 
method to mitigate this jamming is not to accept air¬ 
planes beyond the point at which the airport performance 
decreases. This can be realized by rerouting techniques 
and the Ground Delay Program (GDP) [39j operated by 
the Federal Aviation Administration [40l-l43|. Therefore, 
knowing the flux-density relationship in each airport for 
operational reference is vital. 

Note that the jamming presented in this study is a 
completely different phenomenon from the delay prop¬ 
agation phenomena that have been actively studied 
[J, Wa. I40h 42I |. In these studies, a main cause of delay is 
the crew and passenger connection disruptions. This type 
of delay propagates downstream. In this study, conges¬ 
tion in one airport does not simply propagate; it disrupts 
the density in neighboring airports through its outflow 
traffic. This study also provides a reason why a small 
disturbance, for example, due to a bad weather, grows to 
a large delay from a systemic perspective. Furthermore, 
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note that if a congested airport stops accepting aircraft 
by following a strategy such as the GDP, congestion prop¬ 
agates to upstream airports. Thus, the propagation di¬ 
rection is determined by its cause. 

Our results are based on a model that only considers 
fixed physical constraints (surface congestion curve) and 
ignores flight schedules and operational policies. In the 
future, we intend to address congestion problems on avi¬ 
ation networks more comprehensively by integrating all 


these factors, which we believe will provide insight to the 
problem. 
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